classdef berk_ml
	methods (Static)

	% -------------------------------------------------

	% function estimate
	function estimates=estimate(data,config,constr)

		p=data.p;
		y=data.y;
		q=data.q;
		share=data.share;
		sigma_berkson_vec=data.sigma_berkson_vec; 
		mean_berkson_vec=data.mean_berkson_vec; 
		eps_vec=data.eps_vec;
		eps_long=data.eps_long;

		K1=config.K1;
		K2=config.K2;
		K3=config.K3;

		Berkson_estimated=config.Berkson_estimated;
		unit_TT=config.unit_TT;

		p_lb=config.p_lb;
		p_ub=config.p_ub;

		y_lb=config.y_lb;
		y_ub=config.y_ub;

		q_lb=config.q_lb;
		q_ub=config.q_ub;

		computetime=tic;

		mu1=mean_berkson_vec;
		sigma1=sigma_berkson_vec;
		mu2=NaN;
		mu3=NaN;
		sigma2=NaN;
		sigma3=NaN;
		fdist1='Normal';
		fdist2='';
		fdist3='';
		alpha1=NaN;
		alpha2=NaN;
		alpha3=NaN;
		[TT, ~, TT_deriv2, ~] = construct_chebyshev_berkson_Ponly_twodim_inc(p,q,y,K1,K2,K3,p_lb,q_lb,y_lb,p_ub,q_ub,y_ub,NaN,NaN,NaN,NaN,NaN,NaN,mu1,mu2,mu3,sigma1,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3);
		matrix_aggr=eye(size(p,1));
		
		if Berkson_estimated==1
			
			N_hist_eps=size(eps_vec,1);

			p_ext_eps=kron(ones(N_hist_eps,1),p);
			q_ext_eps=kron(ones(N_hist_eps,1),q);
			y_ext_eps=kron(ones(N_hist_eps,1),y);
			
			sigma1_ext_eps=kron(ones(N_hist_eps,1),sigma1);
			
			mu1_ext_eps=eps_long;
			
			[TT, ~, TT_deriv2, ~] = construct_chebyshev_berkson_Ponly_twodim_inc(p_ext_eps,q_ext_eps,y_ext_eps,K1,K2,K3,p_lb,q_lb,y_lb,p_ub,q_ub,y_ub,NaN,NaN,NaN,NaN,NaN,NaN,mu1_ext_eps,mu2,mu3,sigma1_ext_eps,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3);
			
			matrix_aggr=kron(ones(1,N_hist_eps),eye(size(p,1)));
		
		end

		[BB, BB_deriv1, BB_deriv2, BB_deriv3] = construct_chebyshev_berkson_Ponly_twodim_inc(p,q,y,K1,K2,K3,p_lb,q_lb,y_lb,p_ub,q_ub,y_ub,NaN,NaN,NaN,NaN,NaN,NaN,0*mu1,mu2,mu3,0*sigma1,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3);

			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			% NAG
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

		user.mat1=TT_deriv2;

		if Berkson_estimated==1
			user.mat1=matrix_aggr*TT_deriv2;
		end
		
			% 

		N=size(q,1);
		
		if unit_TT==0
			a=[matrix_aggr*TT_deriv2; BB_deriv2; BB]; 
			sizeTT=size(matrix_aggr*TT_deriv2,1);
			lb=[ -9*ones((K1+1)*(K2+1)*(K3+1),1); 0.001*ones(sizeTT,1); 0.001*ones(N,1); 0*ones(N,1)];
			ub=[  9*ones((K1+1)*(K2+1)*(K3+1),1);   999*ones(sizeTT,1);   999*ones(N,1); 1*ones(N,1)];
		end
		if unit_TT==1
			a=[matrix_aggr*TT_deriv2; BB_deriv2; BB; matrix_aggr*TT]; 
			sizeTT=size(matrix_aggr*TT_deriv2,1);
			lb=[ -9*ones((K1+1)*(K2+1)*(K3+1),1); 0.001*ones(sizeTT,1); 0.001*ones(N,1); 0*ones(N,1); 0*ones(N,1)];
			ub=[  9*ones((K1+1)*(K2+1)*(K3+1),1);   999*ones(sizeTT,1);   999*ones(N,1); 1*ones(N,1); 1*ones(N,1)];
		end

		if Berkson_estimated==0
		
			perc_q=NaN*q;
			parfor ii=1:N
				perc_q(ii)=mean(q<q(ii)+p(ii));
			end
			hat=(TT'*TT)^(-1)*TT'*perc_q;

			theta0=hat*1.01;
		
		end
		if Berkson_estimated==1
		
			perc_q=NaN*q;
			parfor ii=1:size(perc_q,1)
				perc_q(ii)=mean(q<q(ii)+p(ii));
			end
			perc_q=kron(ones(N_hist_eps,1),perc_q);

			hat=(TT'*TT)^(-1)*TT'*perc_q;

			theta0=hat*1.01;
		
		end
		
		if ~constr	

			if unit_TT==0
				N_constr=sizeTT + 2*N;
			end
			if unit_TT==1
				N_constr=sizeTT + 2*N + sizeTT;
			end

			[theta_MLE_unconstr,success_unconstr,fval_MLE_unconstr]=run_NAG(N_constr,lb,ub,a,theta0,user);
			estimates.theta_hat=theta_MLE_unconstr;
			estimates.success=success_unconstr;
			estimates.fval=fval_MLE_unconstr;
	
		end

		if constr	

			sel_impose=(p>=0.2 & p<=0.36 & q<percentiles(q,90) & q>percentiles(q,10) & y>=log(20000) & y<=log(90000));
			N_sel_impose=sum(sel_impose,1);

			a=[a; -BB_deriv1(sel_impose,:) - (share(sel_impose,:)*ones(1,size(BB_deriv3(sel_impose,:),2))).* BB_deriv3(sel_impose,:)];

			lb=[lb; -999*ones(N_sel_impose,1)];
			ub=[ub; 0*ones(N_sel_impose,1)];

			if unit_TT==0
				N_constr=sizeTT + 2*N + N_sel_impose;
				assert(N_constr+size(theta0,1)==size(lb,1))
			end
			if unit_TT==1
				N_constr=sizeTT + 2*N + N_sel_impose + sizeTT;
				assert(N_constr+size(theta0,1)==size(lb,1))
			end

			[theta_MLE_constr,success_constr,fval_MLE_constr]=run_NAG(N_constr,lb,ub,a,theta0,user);
			estimates.theta_hat=theta_MLE_constr;
			estimates.success=success_constr;
			estimates.fval=fval_MLE_constr;
			
		end

		estimates.timer=toc(computetime);
		estimates.config=config;

	end % end of berk_ml.estimate
	% -------------------------------------------------
	
	function estimates=estimate_noY(data,config,sigma_berkson)

		p=data.p;
		q=data.q;
		K1=config.K1;
		K2=config.K2;
		Berkson_estimated=config.Berkson_estimated;

		p_lb=config.p_lb;
		p_ub=config.p_ub;

		q_lb=config.q_lb;
		q_ub=config.q_ub;
		
		mu1=0;
		sigma1=sigma_berkson;
		fdist1='Normal';
		alpha1=NaN;
		[TT, ~, TT_deriv2] = construct_chebyshev_berkson_Ponly_twodim(p,q,K1,K2,p_lb,q_lb,p_ub,q_ub,NaN,NaN,mu1,sigma1,fdist1,alpha1);

		[BB, ~, BB_deriv2] = construct_chebyshev_berkson_Ponly_twodim(p,q,K1,K2,p_lb,q_lb,p_ub,q_ub,NaN,NaN,mu1,0*sigma1,fdist1,alpha1);

			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			% NAG
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

		user.mat1=TT_deriv2;

		if Berkson_estimated>=1
			user.mat1=matrix_aggr*TT_deriv2;
		end

			% 

		N=size(q,1);
		a=[TT_deriv2; BB_deriv2; BB]; 
		lb=[ -9*ones((K1+1)*(K2+1),1); 0.001*ones(N,1); 0.001*ones(N,1); 0*ones(N,1)];
		ub=[  9*ones((K1+1)*(K2+1),1);   999*ones(N,1);   999*ones(N,1); 1*ones(N,1)];
		
		perc_q=NaN*q;
		for ii=1:N
			perc_q(ii)=mean(q<q(ii)+p(ii));
		end
		hat=TT\perc_q;
		theta0=hat*1.01;

		N_constr=3*N;
		[theta_MLE_unconstr,~,~]=run_NAG(N_constr,lb,ub,a,theta0,user);

		estimates.theta_hat=theta_MLE_unconstr;
		estimates.TT=TT;

	end % end of berk_ml.estimate_noY
	% -------------------------------------------------


	function [data,config]=load_data(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_berkson,berkson_off,Berkson_estimated,unit_TT,input_directory,output_directory)
		
		input=load([input_directory 'data_BHP2_trim.mat']);
		p=input.log_p;
		y=input.log_y;
		q=input.log_q;
        distance_oil1000=input.distance_oil1000;
		state_fips=input.state_fips;
		w=(distance_oil1000-min(distance_oil1000))/(max(distance_oil1000)-min(distance_oil1000));
		config.theta_true=NaN;

			% TRIM DATA 
		 
		if trim_perc_q>0.001
			perc_q_low=percentiles(q,trim_perc_q);
			perc_q_high=percentiles(q,100-trim_perc_q);
			sel=(q>=perc_q_low & q<=perc_q_high);
			p=p(sel);
			y=y(sel);
			q=q(sel);
			w=w(sel);
			state_fips=state_fips(sel);
		end

		if min(isnan(y))==0    

			sel=(y>=log(y_min) & y<=log(y_max));

			p=p(sel);
			y=y(sel);
			q=q(sel);
			w=w(sel);
			state_fips=state_fips(sel);

		end
		
			% LOAD BERKSON DATA 

		file_sd=load([input_directory 'berkson_sd_state.mat']);
		tt=[file_sd.state_fips file_sd.sd_uhat];

		sd_berkson=NaN*p;
		mean_berkson=0*p;

		if strcmp(berkson_txt,'state_hist0')==1
			parfor ii=1:size(p,1)
				temp=tt(tt.state_fips==state_fips(ii),2); 
				sd_berkson(ii)=temp{1,1}; 
			end
		end
		if strcmp(berkson_txt,'base_hist0')==1
			sd_berkson=0.033*ones(size(p,1),1);
		end

		if berkson_off==1
			sd_berkson=0*sd_berkson;
		end

		sd_berkson=factor_berkson*sd_berkson;
		eps_vec=NaN;
		eps_long=NaN;

		if Berkson_estimated==1
			
			eps_vec=csvread([output_directory 'station_data_eps_vec_N1000.csv']);
			N_hist_eps=size(eps_vec,1);
			input_mat=csvread([output_directory 'station_data_eps_vec_state.csv']);
			eps_long=NaN*ones(size(p,1)*N_hist_eps,1);
			ii_long=kron(ones(N_hist_eps,1),(1:size(p,1))');
			for ii=1:size(p,1) 
				vec_sel=(input_mat(:,1)==state_fips(ii));
				aux_dd=(ii_long==ii);
				eps_long(aux_dd,1)=input_mat(vec_sel,2); 
			end

			factor_bw=1;
			h_s=factor_bw*1.06*std(eps_vec)*N_hist_eps^(-1/5);
			sd_berkson=h_s*ones(size(p,1),1);

			if berkson_off==1
				sd_berkson=0*sd_berkson;
				eps_vec=0*eps_vec;
				eps_long=0*eps_long;
			end
		end

		share=(p.*q)./y;

		p_lb=min(p)-range_extend;
		p_ub=max(p)+range_extend;
		y_lb=min(y);
		y_ub=max(y);

		q_lb=min(q);
		q_ub=max(q);

		q(q<q_lb)=q_lb;
		q(q>q_ub)=q_ub;

		y_interest_univ=log([72500; 57500; 42500]);
		p_interest_univ=linspace(percentiles(p,5), percentiles(p,95),50)';

		y_interest_univ_DWL=y_interest_univ;
		p_interest_univ_DWL=p_interest_univ; 

		data.p=p;
		data.y=y;
		data.q=q;
		data.share=share;
		data.w=w;
		data.sigma_berkson_vec=sd_berkson; 
	    data.mean_berkson_vec=mean_berkson; 
		data.eps_vec=eps_vec;
		data.eps_long=eps_long;

		config.p_lb=p_lb;
		config.p_ub=p_ub;
		config.y_lb=y_lb;
		config.y_ub=y_ub;

		config.q_lb=q_lb;
		config.q_ub=q_ub;

		config.q0=mean(q,1);

		config.K1=K1;
		config.K2=K2;
		config.K3=K3;

		config.mat_type=mat_type;
		config.N=size(q,1);

		config.p_interest_univ=p_interest_univ;
		config.y_interest_univ=y_interest_univ;

		config.p_interest_univ_DWL=p_interest_univ_DWL;
		config.y_interest_univ_DWL=y_interest_univ_DWL;

		config.Berkson_estimated=Berkson_estimated;
		config.unit_TT=unit_TT;

	end % end of load_data
	% -------------------------------------------------

	function [data_bootstrap,config]=load_data_bootstrap(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_berkson,berkson_off,Berkson_estimated,unit_TT,input_directory,output_directory)

		[data,config]=berk_ml.load_data(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_berkson,berkson_off,Berkson_estimated,unit_TT,input_directory,output_directory);
		
		p=data.p;
		y=data.y;
		q=data.q;
		share=data.share;
		w=data.w;
		sigma_berkson_vec=data.sigma_berkson_vec;
		mean_berkson_vec=data.mean_berkson_vec;

			% DRAW DATA
		
		N=size(q,1);
		[p_bootstrap, y_bootstrap, q_bootstrap, share_bootstrap, w_bootstrap, sigma_berkson_vec_bootstrap, mean_berkson_vec_bootstrap] = resample_data_flex(N,'with',p, y, q, share, w,sigma_berkson_vec, mean_berkson_vec);
		
		data_bootstrap.p=p_bootstrap;
		data_bootstrap.y=y_bootstrap;
		data_bootstrap.q=q_bootstrap;
		data_bootstrap.share=share_bootstrap;
		data_bootstrap.w=w_bootstrap;
		data_bootstrap.sigma_berkson_vec=sigma_berkson_vec_bootstrap;
	    data_bootstrap.mean_berkson_vec=mean_berkson_vec_bootstrap;
		data_bootstrap.eps_vec=NaN;
		data_bootstrap.eps_long=NaN;

	end % end of load_data_bootstrap
	% -------------------------------------------------

	function compute_DWL(data,est1,est1_legend,est2,est2_legend,est3,est3_legend,est4,est4_legend,ref,output_folder,show_txt,addtxt,show_DWL_flag,flag_median)

		q=data.q;
		mu1=0;
		mu2=NaN;
		mu3=NaN;
		sigma2=NaN;
		sigma3=NaN;
		fdist1='Normal';
		fdist2='';
		fdist3='';
		alpha1=NaN;
		alpha2=NaN;
		alpha3=NaN;
		config=est1.config;

		q_lb=config.q_lb;
		q_ub=config.q_ub;
	
		y_interest_univ=config.y_interest_univ;
		p_interest_univ=config.p_interest_univ;

		q0=mean(q,1);   
		p_interest=kron(p_interest_univ,ones(size(y_interest_univ,1),1));
		y_interest=kron(ones(size(p_interest_univ,1),1),y_interest_univ);
		tau_vec=[0.25; 0.50; 0.75]; 

		q_hat_est1=berk_ml.invert_Ginv(est1,q0,tau_vec,p_interest,y_interest);
		q_hat_est2=berk_ml.invert_Ginv(est2,q0,tau_vec,p_interest,y_interest);
		q_hat_est3=berk_ml.invert_Ginv(est3,q0,tau_vec,p_interest,y_interest);
		q_hat_est4=berk_ml.invert_Ginv(est4,q0,tau_vec,p_interest,y_interest);

		config=est1.config;
		K1=config.K1;
		K2=config.K2;
		K3=config.K3;
		p_lb=config.p_lb;
		p_ub=config.p_ub;
		y_lb=config.y_lb;
		y_ub=config.y_ub;
		
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			% COMPUTE DWL
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
		if show_DWL_flag==1

			fname=[output_folder 'DWL_' ref '.txt'];
			fileID=fopen(fname,'w');
			fprintf(fileID,'tau;  inc; DWL_pertax; DWL_perinc; ref\n');
			fclose(fileID);

			show_DWL(p_interest_univ,y_interest_univ,p_interest,y_interest,tau_vec,K1,K2,K3,p_lb,p_ub,q_lb,q_ub,y_lb,y_ub,mu1,mu2,mu3,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3,est1,q_hat_est1,show_txt,[est1_legend '; ' addtxt],fname,flag_median);
			show_DWL(p_interest_univ,y_interest_univ,p_interest,y_interest,tau_vec,K1,K2,K3,p_lb,p_ub,q_lb,q_ub,y_lb,y_ub,mu1,mu2,mu3,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3,est2,q_hat_est2,show_txt,[est2_legend '; ' addtxt],fname,flag_median);
			show_DWL(p_interest_univ,y_interest_univ,p_interest,y_interest,tau_vec,K1,K2,K3,p_lb,p_ub,q_lb,q_ub,y_lb,y_ub,mu1,mu2,mu3,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3,est3,q_hat_est3,show_txt,[est3_legend '; ' addtxt],fname,flag_median);
			show_DWL(p_interest_univ,y_interest_univ,p_interest,y_interest,tau_vec,K1,K2,K3,p_lb,p_ub,q_lb,q_ub,y_lb,y_ub,mu1,mu2,mu3,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3,est4,q_hat_est4,show_txt,[est4_legend '; ' addtxt],fname,flag_median);

		end

		ghat_MLE_tau50_unconstr=q_hat_est1(:,2);
		ghat_MLE_tau50_constr=q_hat_est2(:,2);
		ghat_MLE_tau50_unconstr_noBerk=q_hat_est3(:,2);
		ghat_MLE_tau50_constr_noBerk=q_hat_est4(:,2);

		ghat_MLE_tau25_unconstr=q_hat_est1(:,1);
		ghat_MLE_tau25_constr=q_hat_est2(:,1);
		ghat_MLE_tau25_unconstr_noBerk=q_hat_est3(:,1);
		ghat_MLE_tau25_constr_noBerk=q_hat_est4(:,1);

		ghat_MLE_tau75_unconstr=q_hat_est1(:,3);
		ghat_MLE_tau75_constr=q_hat_est2(:,3);
		ghat_MLE_tau75_unconstr_noBerk=q_hat_est3(:,3);
		ghat_MLE_tau75_constr_noBerk=q_hat_est4(:,3);

		save([output_folder 'data_' ref '.mat'],'p_interest*','y_interest*','ghat*');
		
	end % end of compute_DWL
	%--------------------------------------------------
	
	% function invert_Ginv
	function q_hat=invert_Ginv(est1,q0,tau_vec,p_interest,y_interest)

		config=est1.config;
		theta_hat=est1.theta_hat;
		q_hat=NaN*p_interest*tau_vec';

		options = optimoptions('fsolve','Display','off');
					
		if max(isnan(theta_hat))==0
		
			if length(tau_vec)>1
				for pp=1:length(p_interest)
					p_interest_sel=p_interest(pp);
					y_interest_sel=y_interest(pp);
					for tt=1:length(tau_vec) 
						tau=tau_vec(tt);
						[q_hat(pp,tt),~,~]=fsolve(@(q) find_q_with_inc(p_interest_sel,q,y_interest_sel,tau,theta_hat,config),q0,options);
					end
				end
			end
			if length(tau_vec)==1
				tau=tau_vec(1);
				parfor pp=1:length(p_interest) 
					p_interest_sel=p_interest(pp);
					y_interest_sel=y_interest(pp);
					[q_hat(pp,1),~,~]=fsolve(@(q) find_q_with_inc(p_interest_sel,q,y_interest_sel,tau,theta_hat,config),q0,options);
				end
			end

		end

	end % end of invert_Ginv
	%--------------------------------------------------

	% function invert_Ginv_noY
	function q_hat=invert_Ginv_noY(est1,q0,tau_vec,p_interest)

		config=est1.config;
		theta_hat=est1.theta_hat;
		q_hat=NaN*p_interest*tau_vec';

		options = optimoptions('fsolve','Display','off');
					
		if max(isnan(theta_hat))==0
		
			if length(tau_vec)>1
				for pp=1:length(p_interest)
					p_interest_sel=p_interest(pp);
					parfor tt=1:length(tau_vec) 
						tau=tau_vec(tt);
						[q_hat(pp,tt),~,~]=fsolve(@(q) find_q_noY(p_interest_sel,q,tau,theta_hat,config),q0,options);
					end
				end
			end
			if length(tau_vec)==1
                tau=tau_vec(1);
				parfor pp=1:length(p_interest) 
					p_interest_sel=p_interest(pp);
					[q_hat(pp,1),~,~]=fsolve(@(q) find_q_noY(p_interest_sel,q,tau,theta_hat,config),q0,options);
				end
			end

		end

	end % end of invert_Ginv_noY
	%--------------------------------------------------

	function main_compute(berkson_ref,Berkson_estimated,K1,K2,K3,ff,unit_TT,input_directory,output_directory)

		berkson_txt=[berkson_ref '_hist' num2str(Berkson_estimated)];
		factor_vec=[1; 1.2; 0.8; 0.5];
		range_extend=0.033; 
		mat_type='chebyshev'; 
		trim_perc_q=1;
		y_min=1;
		y_max=999000;

		fprintf('--------------------------------------------------------------\n');
		fprintf(['berkson_txt=%5s, K1=%3.0f, K2=%3.0f, K3=%3.0f, ff=%1.0f [' datestr(clock) ']\n'],berkson_txt,K1,K2,K3,ff);
		fprintf('--------------------------------------------------------------\n');

		[data,config]=berk_ml.load_data(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_vec(ff),0,Berkson_estimated,unit_TT,input_directory,output_directory);
		[data_noBerk,config_noBerk]=berk_ml.load_data(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,berkson_txt,factor_vec(ff),1,Berkson_estimated,unit_TT,input_directory,output_directory);

		estimates_unconstr_noBerk=berk_ml.estimate(data_noBerk,config_noBerk,0);
		estimates_unconstr=berk_ml.estimate(data,config,0);
		estimates_constr=berk_ml.estimate(data,config,1);
		estimates_constr_noBerk=berk_ml.estimate(data_noBerk,config_noBerk,1);

		estimates_filename=['workfile_' berkson_txt '_ff' num2str(ff) '_K1_' num2str(K1) '_K2_' num2str(K2) '_K3_' num2str(K3) '_' berkson_txt '.mat'];
		filename=[output_directory estimates_filename];	
		save(filename);

	end % end of main_compute
	%--------------------------------------------------
							
	function main_output(berkson_ref,Berkson_estimated,K1,K2,K3,ff,output_directory,show_figs,show_combined,show_income,show_DWL_flag,flag_median)

		berkson_txt=[berkson_ref '_hist' num2str(Berkson_estimated)];
		addtxt=sprintf('K1=%3.0f, K2=%3.0f, K3=%3.0f, txt=%5s',K1,K2,K3,berkson_txt);
		file_ref=[berkson_txt '_ff' num2str(ff) '_allY_p' num2str(K1) '_q' num2str(K2) '_y' num2str(K3)];
		ref_txt=['ff' num2str(ff)];

		filename=[output_directory 'workfile_' berkson_txt '_ff' num2str(ff) '_K1_' num2str(K1) '_K2_' num2str(K2) '_K3_' num2str(K3) '_' berkson_txt '.mat'];	
		
		S=load(filename);
		data=S.data;
		estimates_unconstr=S.estimates_unconstr;
		estimates_constr=S.estimates_constr;
		estimates_unconstr_noBerk=S.estimates_unconstr_noBerk;
		estimates_constr_noBerk=S.estimates_constr_noBerk;
	
		berk_ml.compute_DWL(data, ...
					estimates_unconstr,'unconstr', ...
					estimates_constr,'constr', ...
					estimates_unconstr_noBerk,'no Berkson, unconstr', ...
					estimates_constr_noBerk,'no Berkson, constr', ...
					[berkson_txt '_ff' num2str(ff) '_allY_p' num2str(K1) '_q' num2str(K2) '_y' num2str(K3)], ...
					output_directory,1,addtxt,show_DWL_flag,flag_median); 

		if show_figs==1
			prepare_graphs(file_ref,ref_txt,output_directory,show_combined,show_income);
		end

	end % end of main_output
	%--------------------------------------------------
							
	%
	end % end of methods
	% -------------------------------------------------

end % end of berk_ml class
% -------------------------------------------------

	%
